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Results to Date 

We have now received two Landsat*4 TM tapes (7 bands) within our study area in 
the southern Sierra Nevada (path 41, row 3b, 10 December 1962, and path 42, row 34, 
13 January 1963). During the reporting period we have worked on registration of the 
TM data to digital topographic data, on oomparlion of TM, MSS and NOAA meteorologi- 
cal satellite data for snowcover mapping, and on radiative transfer models for atmos- 
pheric correction. 

These investigations are described in more detail below. 



^^Account ledger! for June are incomplete. 
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PresentatiMU by Princi|Ml Invasticator 

1. April 20, 1963, Snow Rtfltctanct from tha Thomaiic Mapper, Institute of Remote 
Sensing, Academia Sinica, Beijing, China. 

2. i^rii 20, 1963, Slope, Horvean, and Draimtge ^isin fin/ormaHon from JXgiial Ter- 
rain Models, Institute (tf Remote Sensing, Academia Sinica, Beijing, China. 

3. April 2Z, 1963, Remote Sensing and Snow Hydrology, Institute of Remote Sensing, 
Academia Sinica, Beying. China. 

4. April 22, 1963, Image Processing Software at Uniuersity of CMifomia, Santa Bar- 
bara, Institute of Remote Sensing, Acadeicia Sinica, Beijing, China. 

5. April 26, 1963, Emisswity and Infrared Brightness Temperature of Snow, Institute 
of Remote Sensing. Academia Sinica, Beying, China. 

6. May 1, 1963, Snow Hydrology Research and Remote Sensing, Hen Shan Glaciologi- 
cal Station, UrOmqi River ^in. China. 

7. May 4, 1963, Snow Hydrology Research and Remote Sensing, Xinjiang Institute of 
Geography, DrOmc^, China. 

6. May 5, 1963, Punjdamsntal Ra/tomstric Properties, Xityiang Institute of Geogra* 
phy, UrOmqi, China. 

9. May 5, 1^, Emisswity of Snow, Xinjiang Institute of Geography, UrOmqi, China. 

10. May 7, 1963, Remote Sensing and the Snow Surface Radiation Budget, Institute of 
Glaciology and Cryopedolc^, Academia Sinica, Lanzhou, China. 

11. May 7, 1983, Bmissiuity of Snow, Institute of Glaciology and Cryopedology, 
Academia Sinica. Lanzhou. China. 

12. May 6, 1963, Remote Sensir^ tn the U.S., Gansu Province Institute of Remote 
Sensing, Lanzhou, China. 

13. June 6. 1963, Remote Sensing of the Snow Surface Radiatwn Budget, California 
Space Institute, Scripps Institute of Oceanography, University of California, San 
Diego. 


Detailed Descriptions of Investigations 

1. EFTECT OF SPATIAL RESCMAJnm ON THE ACCURACY OF MAPPING SNOV COVERED 
AREA 

1.1, SJ^mopais 

Satellite data from three sensors. Landsat-4 TM and MSS, and NOAA-7 AVHRR, will 
be used to map snow covered area (SCA) over selected basins to determine the error 
associated with using coarser spatial resolution. This information is needed to quantify 
the information lost when high resolution data are not available and to determine the 
minimum resolution requirements for basins of different sizes. The effect of a discon- 
tinuous snowcover within a basin will also be evaluated. This effort will improve our 
ability to make quantitative comparisons of surface measurements made with different 
satellite instruments 

Satellite remote sensing has become increasingly important to hydrologists 
because it provides information on the spatial distribution of parameters of hydrologic 
importance. In snow and ice studies, remote sensing has been used to improve &e 
monitoring of existing conditions and has been incorporated into several runoff fore- 
casting and management systems. Rango [1979] presents a detailed review of 
research in this area throxigh 1976, and the Snow and Ice, chapter of the Sth Pecora 
Symposium [Deutsch et al., 1961] has descriptions of several applications currently 
implemented or being tested. 



-3- 


The moat common (^eroUonal use of remote sensing in snow studies is to monitor 
the areal extent of the snowcover [Bamas and Bbwlay. 1974; fbstar and Jtango, 1975; 
JlfcOtnnis at al, 1975; Rango and Ittan, 1976; Range and Salomanson, 1976X Thrae 
efforts led other investigators to eoDclude that satellite derived estimates of snow 
covered area [SCA] could be used to improve forecasting of snowmelt runoff [Hartna- 
ford, 1977; Range at oi., 1977; Range, 1976; Range at al„ 1979; Brown at al., I960]. 
Uartinac [1975] and Range and Martinae [1979] carried this a step further in their 
model by including satellite derived measurements of SCA as an index to melt. Satel- 
lites can also be viewed as space-borne radiometers. Techniques have been developed 
to use these data to measure or estimate snow surface characteristics that are neces- 
sary for calculation of the surface radiation budget [Dowiar, 1960a. 1961, 1962, 1963; 
Dozier at al.. 1960, 1961; Dozier and Frew, 1961; Dozier and Warren, 1962; Pramjtton 
and Marks, 1960; Marks, 1962]. While these efforts go beyond mapping SCA. any use of 
satellite data for snow l^drology begins with the mapping of the snow cover distribu- 
tion over the surface. To determine the spatial accuracy of any type of satellite data, 
we begin with the estimate of SCA. 

The efforts cited above have all been influenced by the availability, cost, and 
characteristics of data from satellite instruments that are currently operational. 
Operational s|>plications have usually been limited to tl^ measurement of SCA from 
uncorrected photographic prodiu:ts using manual photo interpretation techniques. 
IMgital snow mapping techniques, while faster and less subjective, have been limited to 
specific experiments because of processing difficulties and the high cost of satellite 
data in digital form. For these reasons no effort has been made to assess the relative 
spatial accuracy of utilizing satellite data in different forms or from different instru- 
ments to estimate SCA. The p\irp<»e of this investigation is to conduct an experiment 
to determine the advantages and disadvantages of using data from different satellite 
instruments. From this we will evaluate the effect of basin size and shape, and discon- 
tinuous snowcover on the accuracy of our estimate of SCA for several different types of 
satellite data. 

1.2. SateUiteDaU 

Environmental satellite data are available in a variety of forms that are useful to 
snow hydrologists. For our investigations on the effect of resolution we will use data 
from the Landsat Thematic Mapper (TM) and Multispectral Spectral Scanner (MSS), and 
from the NOAA Advanced Very High Resolution Radiometer (AVHRR). These are mul- 
tichannel instruments that are frequently used to map SCA o^'er basins from 10A:m^ to 
10,000 Table 1 presents a comparison of these instruments. 
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Table 1 




OiaraetariBtlea of TV. HS^ and AVHRR Shiite Dala 



(spectral, spatial, temporal, and economic) 


J/mdsat-4 Thunatic JHapptr (TM) 


waualangths 

nomtnoi 

overpass 

digital data 

band 

Hi 


vixal siaa 

vntarual 

cost 

n 

.452 

- 

.516 

26.5 m 

16 days 

S2500 


.529 

- 

.610 



WM 

.624 

- 

.693 





.776 

- 

.905 





1.568 

• 

1.784 




H 

2.097 

- 

2.347 




HI 

10.422 

_ 

11.661 

114 m 



Landsat-4 MvUisjmctral Scannar (MSS) 


.497 

- 

.607 

57 m 

16 days 

S650 


.603 

- 

.697 




.704 

- 

.814 







- 

-1.036 




NOAA-7 Advancad Vary High Rasolution Radiometar (AVHRR) 

1 

.56 

- 

.72 

1.1 km 

12 hrs. 

S75 

2 

.71 

- 

.98 




3 

3.53 

- 

3.94 




4 

10.32 

- 

11.38 




5 

11.45 

- 

12.42 




Digital data costs are 

approximate for the TM as they are not yet available 

to the public. For this investigation we have the TM data in hand. Digital 
data costs for the AVHRR reflect our costs to travel to Scripps SSOF, use of 

! the facility, and return to UCSB. 

It is assumed that at least four data sets 

are acquired per trip. 






We have acquired TM and MSS data for the southern Sierra Nevada from NASA God- 
dard for December 10, 1082. We have received an additional TM data set for the central 
Sierra Nevada for January 13. 1983. We will acquire AVHRR data from the Scripps Satel- 
lite Oceanographic Facility (SSOF) for the same overpass dates, and we hope to acquire 
one additional MSS data set to correspond to the January TM overpass. 
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1.3. Ibe Eip«rliii«it 

By acquiring data from three different environmental satellites over the same area 
on the same date, we have the uncommon opportunity to compare the spatial accuracy 
associated with using each type to estimate SCA. We plan to answer the following fun* 
damentat questions: 

« For basins of several different sizes and shapes, how much difference is there in 
our estimate of SCA when using TM, MSS, and AVHRR data? 

* Under vdtat conditions of basin size, and snow cover distribution will coarse resolu* 
tion AVHRR data provide adequate information on SCA? 

» How can occasional inputs from TM or MSS data be used to improve estimates of 
SCA that are regularly made with AVHRR data? 

These questions are important because while high resolution TM data are undoubtedly 
more accurate, they are expensive, difficult to acquire with reasonable regularity, and 
be3rond the digital processing capability of most researchers. MSS data are the most 
commonly used satellite data, but they also have poor temporal resolution, and are 
expensive. AVHRR data are inexpensive, and the satellite has frequent overpasses, but 
the coarse spatial resolution makes the data difficult to use in small basins, or under 
conditions of sparse or patchy snowcovers. 

The research is divided logically into two parts. The first step is to calibrate and 
register the three satellite data sets to a standard reference. This is by far the most 
difficult part of the work, and will take the most time and effort. The development of 
the QDIPS image processing system at the Computer Systems Laboratory, University of 
California, Santa Barbara, makes this registration possible. The second involves the 
selection of test basins of varying sizes and snowcover distributions for the mapping of 
SCA using each type of satellite data. Each of these is described in detail below. 

1.3.1. CclUiratmg and Ragataring SataUiia Data Sats to a Refarance Orid 
Before the satellite data can be analyzed, the following processing is required; 

» Radiometric correction; Pixel values must be converted to the appropriate physi- 
cal units (siurface exitance or brightness temperature). 

• Geometric correction: The satellite imagery must be geometrically transformed so 
that it will overlay our digital terrain database. 

The flow of processing is as follows: 

1.3.2. RtuHomatric correction 

TM, MSS, and AVHRR shortwave data are converted to surface exitances using a 
simple linear transform whose multiplicative and additive terms are assumed to be 
constant [Lourifson at at., 1979; AVigef, 1960; Clark and Dasgupta, 1963]. Since these 
coefficients may change over an image frame time, the radiometric correction must be 
performed before any geometric corrections alter the original time-space relationships 
within the image. 

1.3.3. Oiomtfric corrscfion 

Our terrain database is constructed from USGS digital elevation data. Ihe data 
are supplied in geodetic coordinates (latitude, longitude). During the process of 
extracting a region of interest, we transform the data to the Universal TVansverse Mer- 
cator (UTM) projection. We have adopted the UTM projection as our geographic refer- 
ence since it maintains almost constant scale over areas of the sizes in which we are 
interested, and is readily indexed in common units of measurement (meters). Accu- 
rate. computationally efficient formulas are available to convert between UTM and geo- 
detic (latitude and longitude) coordinates [Dozier, 1980b]. The spatial resolution of the 
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individual terrain data sets ("terrain grids") will be resampled to 26.5 meters, 
corresponding to the highest satellite resolution (TM) to be used in this investigation. 

Our problem is therefore to transform the TM, MSS, and AVHRR data into the UTM 
coordinate system of the terrain grids, at the selected spatial resolution. For TM and 
MSS data, the standard digital pr^ucts 8 t 4 >plied by NASA in the Space Oblique Merca* 
tor (SOM) projection. The SOM projection has been mathematically defined [^tyder, 
1981], but computer code implementing the conversions between SOM and geodetic 
coordinates is not currently available. We will therefore have to develop programs for 
mapping from SOM into UTM. Scale changes will be accommodated by resanq>ling. Any 
misregistration will be corrected by polynomial fitting of ground control points (GCPs), 
Experience of other users of TM and MSS data leads us to believe that any such residual 
misregistrations are largely translational in nature, and therefore easily corrected by a 
few GCPs [Dauid PUts, NASA Johnson Space Flight Center, personal communication]. 

AVHRR data as obtained from the Satellite Oceanographic Facility at the Scripps 
Institute of Oceanography have had no geometric corrections performed. Ibe general 
problem of georeferencing satellite imagery requires precise determination of the 
position (altitude and nadir point) and attitude (yaw, pitch, and roll) of the satellite, 
together with a model of the spatial resolution of the sensor system. We are currently 
Implementing systems to predict satellite positions using a posteriori, orbital element 
data, and to build a satellite attitude vector firom telemetered yaw, pitch, and roll data. 
These models are tuned by least-squares adjustment of their coefficients using control 
points with known geodetic locations. This procedure will be used to determine the 
mapping between AVHRR image coordinates and UTM coordinates, using control points 
selected from pseudo-image displays of the terrain grid. 

The end product of the geometric correction process will be TM, MSS, and AVHRR 
data sets registered to digital terrain models and resampled to 26.5 meter resolution. 
Since this resampling represents an interpolation of up to 1 Ox (in the AVHRR case), we 
will investigate whether high-order resampling techniques affect SCA determination in 
any systematic fashion (we suspect that the effect will be negligible). 

1.3.4. Mapping Snowcovered Area 

Several basins that are partially snowcovered and vary in size from 25 to 2500 
will be selected. To find partially snowcovered basins that are small we will select areas 
near the Owens valley where small .?asins have a larger elevation change. Estimates of 
SCA for each basin will be made using each satellite data t^e. Basins will be digitally 
isolated using a procedure developed by Marks at al., [1BB3J and mapping of SCA will be 
done with existing software in the QDIPS system. In each case, the TM data will be used 
eis the reference standard to determine the "best" estimate of SCA as well as the spa- 
tial contiguity of the snowcover. Results from the MSS and AVHRR data will then be 
compared to this standard. 

In addition to the digital mapping of SCA for each test basia James Foster, an 
experienced photo interpreter from NASA Goddard Space Flight Center will manually 
map SCA for each satellite data type. We will use this information to evaluate the rela- 
tive accuracy of manual over automated snow mapping techniques for satellite data of 
different spatial resduUons. 

1.4. Sig nlflc a n ce of this Res e ar ch 

The immediate importance would be in the improved utility of satellite data for 
mapping SCA. A reduced dependence on more e]q>ensive high resolution data from the 
MSS or TM sensors would save costs, and implementation of automated digital tech- 
niques would improve the efficiency of state and federal water management programs. 
In the long run, the improved understanding of satellite data and their potentials and 
limitations, should provide broad benefits to both hydrology and remote sensing. 
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2. A1V0SPHER1CRM)UT1VET1MNSFER MODEL ACCOUNTING FOR BOni THE ZENITH AND 
AaHUTH VARlATiaNS IN THE RADIAITVE HEU) 

2.1. Synopsis 

In remote sensing, the signature obtained by the satellite sensor is related to radi- 
ance, which is a function of elevation, location, and viewing angle, in addition to the 
astronomic and atmospheric conditions. 

In order to infer the radiometric properties of the medium underlying the atmo- 
sphere from remotely sensed data, atmospheric correction of the satellite data is 
necessary. In remote sensing for oceanography, quick atmospheric correction algo- 
rithms have been developed [6brion, 1976; Gbnton and Clark, 1960] and successfully 
applied to the quantitative assessment of chlorophyll concentration in Southern Cali- 
fornia [,9m.ifA and IKUson, 1961]. In these algorithms, the known spectral optical pro- 
perties of water are used to determine the aerosol path radiance at wavelength 670 
nm, where the subsurface upwelling radiance is negligible, and then to estimate the 
aerosol path radiance at other wavelengths. However, in other wavelengths, or for 
other types of surfaces besides water, the appllcaUon of these algorithms may be lim- 
ited, because of the different optical properties. Therefore, some general atmospheric 
correction algorithms are required. 

In our project, a different approach, which invokes an atmospheric radiative 
transfer model accounting for both the zenith and azimuth variation in the radiative 
field, is taken. The model deals with a plane parallel structured atmosphere, which is 
composed of different layers with each assumed to be homogeneous in composition and 
to have linear-in-r temperature profile. Moreover, the multiple scattering between 
such an atmosphere and the underlying medium (say, water, snow, or soil), which may 
have both the zenith and azimuth variation in radiative behavior, is taken into con- 
sideration. 

This approach is an extension of an azimuthaliy integrated atmospheric radiative 
transfer model developed by H^tscombs [1961], which is based on the results of earlier 
researchers [Qrant and Hunt, 1969] and his own investigations [f^scomAe, 1976a, 
1976b]. The additional parameter on azimuthal variation is handled by adding a new 
azimuth dimension into the computation space. In order to save computation time, 
the Fourier transform is performed over the azimuth domain. Then the resulting 
Fourier coefficients with different order stand for the original vector and work in a 
parallel manner, so that the array sizes for each order of coefficient in the computa- 
tion remain the same as for the azimuthaliy integrated case. The final radiance are 
obtained by inverse Fourier transformation from the Fourier coefficients of thocj quan- 
tities, which are the results derived from the procedures similar to what is used for 
calculating the azimuthaliy integrated radiance. 

The inputs of such a model are some astronomic and atmospheric parameters for 
each layer (or level), and the surface radiative property. Those parameters or proper- 
ties can be obtained either by some model or by measurements. The astronomical 
parameters are earth-sun distance and solar fiux at the top of the atmosphere, and the 
atmospheric parameters Include pressure temperature, cheihical composition of the 
air molecules, and the composition and size of the aerosol, water droplets, and ice cry- 
stals. The outputs of the model are the monochromatic radiance and irradlance of at 
each level. By integration, the total radiance and irradlance can be obtained. 

The potential application of this model in the atmospheric correction of remote 
sensing include several aspects. For example: after calculating the upward radiance at 
the top of the atmosphere for an atmosphere with black surface beneuth it, and then 
by removing this value from the radiance measured at the top of the atmosphere, one 
can get the modified radiance which accounts for the contribution of upward radiance 
at the surface level. 
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Also, this model can mimic the angular radiative field for a plane parallel layered 
atmosphere * earth system, in which the earth surface property is homogeneous from 
location to location. For low frequency spatial variation, this model may also offer good 
approximation. For more complicated cases with high frequency spatial variation, this 
model may approximate the average radiance over an area with certain size. The 
departure from the average may be used to derive more detailed information of the 
surface. 

Besides, this model will be useful to assess the other atmospheric correction algo- 
rithms, in which some aspects of the interaction between radiation and medium may 
be omitted. Because this model considers all the aspects of the interaction between 
photon and medium for a plane parallel structured absorbing and scattering system, if 
we apply the other algorithms on such a system, the departure of the results of the two 
models will give an indication of the performance of the algorithm to be assessed. 

2.2. Review of Radiatitm Model 

Once the atmosphere is subdivided into many plane parallel homogeneous layers 
according to some atmosphere model or some measurement data, the basic inherent 
radiative properties, such as the total optical depth, the single scattering albedo, and 
the phase function can be derived for each layer by Rayleigh scattering. Mie scattering, 
line absorption and continuum absorption theories. 

Based on those basic properties, the reflection and transmission functions, and 
the solar and thermal source function for each layer can be determined, and finally, 
the radiances at desired levels can also be calculated, using the radiative transfer 
theory. 

In the following sections, we will deal with a /T-layer atmosphere plus underlying 
medium model. 

2.2. 1. Radiative Transfer Ejuation 

The radiative transfer equation for a plane parallel medium is given by [Chon- 
drasekhar, 1960] 


T is optical depth, and L{T,fx,^) is the radiance at level r along direction fx, <p (where /x 
is cosine of zenith angle and <p is azimuth angle). The source function is given by 

a i 

- T=-f f L{r,tx',tp')dfx'd(fi' + Q{T.fi.(p) 

w 0 

The last term represents an internal source. The internal source term in this case is 
similar to that for the azimuthally averaged case [ Wiscombe, 1976]. 

Qir.fx,^) = Qt(T.n,^) + 

The thermal source is 

= (l-B)S(f(T)) 

is the Plank function at temperature T, where 7 is a function of optical depth. 
The direct source and specular "pseudo-source” arising from the diffuse - direct split- 
ting of the radiance are 
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Here fio is the cosine of the solar zenith angle, Fq is the solar irradiance incident on the 
top of the atmosphere (normal to the beam), and p, is the directional specular 
reflectivity at the surface beneath the atmosphere. 


2.2.2. IrUaraction Principlt 

One way to attack this integro - dillerential equation is using the "interaction prin- 
ciple” [Ohmf and Hunt, 1969] which, specified to a layer bounded by T/ and Tj across 
which D and phase function P are constant, assumes the discrete form 

= r(T/,T/)L'(T/) + t{rj,r,)L*{rj) + E^(t/,T/) 

= r(Ty.T/)L*(T/) + f (T/.T/)L-(T/) + E-(T;.T/) 


Radieuice L*(t), (tsT/,T/) is a vector of mxn elements on a discrete angular space 
composed of m zenith and n azimuth angles 


L* (f) = 






0</i,< are a set of quadrature points for (0,1). 0s^i< <^n<2rr are 

equally spaced points in the interval 0-2rr. r(T/,T/) and f(T/,T/) are reflection and 
transmission matrices. £*(0,r) are internal source vectors. The r, t and £ for each of 
the homogeneous single layer can be derived by some initialization scheme and the 
doubling method. 


2.2.3. Initialization 

In a recent letter, Wiscombe suggests that the Infinitesimal Generator Initializa* 
tion (IGl) is best [Orant and Hunt, 1969]. It produces the r(Ar), f (At), and S^(Ar) for 
an optically very thin starting layer with optical depth At. After initialization, the 
’■(t/.T/*.|), f (t/,t/^i), and E*(t/,T/ 4 .i) for the entire homogeneous layer bounded by lev- 
els / and I+l can be derived. 

The IGl claims 

t(At) = p*- c 

4TT 

f (A t) -I - I ^ SiSL^^-i p** c 

4TT 

where I.U.C are matrices of size mn xmn. W* is the inverse ot M. / is the identity 
matrix, M and C are diagonal matrices, 

* [iUi d«] 

C - [c< <5y 6«] 

and 

where \ni<m, Kk^n, and Kl^n, and fit's are the cosines of zenith angles 

derived by some quadrature rule (say, Gaussian quadrature rule), and the c^'s are the 
associated weights, and is delta function defined as d^ = l if t=;, otherwise d^=0. 
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The source vectors for dlllerent source need to be handle separately. The IGl gives 
the general relation between source term &nd the source vectors for the initial 
layer 


2*(At) = AT* 


f 


V 


9(t.±Mi.«Ps) 


dr 


Qir,±^.9n) 


The solar source vector for the initial thin layer is 

E| (dr) = — U “ • 
and the thermal source term 


E/(At) = (1-D)^ 

However, the thermal source vectors for the entire homogeneous layer with linear-ln- 
tau temperature profile within it can be derived in a quite simple fashion vrithout 
directly using £^(Ar). 


2.2.4. Adding and Dovbling litthod 

By applying the interaction principle to two adjacent layers, the reflection and 
transmission matrices and the source vectors for the combined layer can be derived if 
the corresponding quantities are known for each of these two layers [f^anf and Hunt, 
1969]. 

Consider two adjacent layers bounded by planes x.y.z. By the interaction princi- 
ple, we have 

Z,*(y) = t{y,x) L*{x) + r(i,y) L~{y) + E*(z.y) 

L~{x) = r(x,y) L*{x) + f (x,y) L~{y) + E‘(x,y) 

L*{x) -t{z,y) L*(y) + r{y,x) L~iz) + E*(y,x) 

i‘(y) = r{y.z) L*{y) ^ t{y,z) L~{z) + E‘(y,s) 

Since x,y,s are entirely arbitrary, we can use the interaction principle and write 

L*{z) - t(z.x) L*{x) + r(x.s) L~{z) + E*(x,s) 

L~{x) = r(x,s) L*{x) + f (x.s) L~{z) + E“(x,*) 

These equations can be obtained by eliminating L*{y) and L~(y) from those before. 
Then the following relations hold 

r(z,x) =r(y,x) + f(x,y) [/ - r(s.y) r(x,y)]‘‘ r(s,y) f (y,x) 

r(x,x) = r(y,s) + f(s,y) [/ -r(x.y)r(*,y)]-‘r(x,y) f(y,z) 

f(z,x) = f(x,y) [/ -r(x,y)r(s,y)]-‘ f(y,x) 

f(x,s) a <(*,y) [/ -r(f .y)r{x,y)]-‘ f{y.f) 

E*^ (x,x) = t{z.y) [/ -r(x.y)r(s,y)]-‘ r(x,y) E*(y,z) + 

s^(y.*) + £*(*.y) 

E'(x,s) a f (x,y) [/ -r(s.y) r(x.y)]-‘ r(s.y) E*(y,x) + 
t{x.y) [I -r{z.y) r(x,y)]-‘ E-(y.s) + E"(x,y) 

These are the formulae for the so called "adding" method. If the two layers are 


• 11 - 


ORIGINAL PAGE tS 
OF POOR QUALITY 


identical, the method is called "doubling." If the initial layer is chosen )$uch that 

At = (T/t, 

where is a integer, and (r/^i -r/) is the optical depth of the I*l th layer in the 
multl'layer system, then the reflection and transmission matrices and source vectors 
for such a homogeneous thick layer can be built up quickly by "doubling" N times: 

fntl “ fn (/ “ T|, r,i) * 

T«ti = + fn (/ - Tn T^)-* T» 

r„=r(i?*flT), fn=f(2"Ar), eiidQ£n*N 

Since neither of the Internal sources, including the solar and thermal sources, are 
constant with optical depth, nor the changes of them with optical depth the same, they 
need to be treated separately. For the solar source, the doubling scheme takes a spe- 
cial form; 

” f|* Tn (•» T„ + £/,n) + Sn 
^e.n + 1 “ ffi (Tn ^^.n ^ *n Ss.n) + 

For initial layer. nsO. ®nd 

For a layer with linear-in-r temperature profile, the doubling method gives the 
flnal resulting thermal source vectors as [ lyiscombs, 1976]: 

Yn±B‘Zn 

where 

B^^B[nr,)] 

BN = B[nrt,^)] 

B' = [{Bs-B^)/{ti,,-t,) 

Zn can be obtained recursively; 

^n*\ ” 9n (^n I'n ~ ^n r*n Tn) {Zy^ — Pn Yf^) 

yn = 2"-'flT 

Yn.l ={tn r« Tn + T« + /) F« 

With the initial values Z^-O, Po=JiAT and Fo“(1“0)At M~W. U is a vector of ones, and 

r=(/-Tnrn)-‘ 

Then the total source vectors for the entire homogeneous layer are the sum of all 
the different types of internal sources. 

+ • • • 


2.2.5. Calculation of the Internal RadUmee 

Having known the reflection and transmission matrices, and the source vectors for 
each layer in the multi-layer system, we can build the internal radiance field in the 
atmosphere by the adding method. Using formulae of the interaction principle, we 
have a set of simultaneous equations: 

A*(T/*i) = f (T/*I.T/) L*{Tj) + r(T/.T/*i) + E*(T/.T/*,) 

L~{ri) = t(t/*,.t/) L*(ri) + 

for (XKK. where K is the total number of layers in the system, and to= 0 and tx is 
total optical depth of the atmosphere. In addition, two sets of boundary conditions 
need to be sa'isfled: 
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L*(to)=L*(0) 

It the top boundary condition. 

I-(Tjr) = Rg i"(Tjr) + * BiTe) + frM 

is the bottom boundary condition. Ro it the surface diffuse reflection )tiatrix, Jr{ja^ is 
the BRDF (bidirectional reflectance-distribution functions) vector for the incident 
beam, and To is the temperature of the surface. 

It is convenient to solve this set of equations in the following recursive manner. 

1. Compute recursively the intermediate partial interactive quantities from top to 
bottom for each layer (or level). 

This is so called forward pass. In this stage, the calculation starts at the top layer, 
then goes down, and for each layer only the contribution made by the internal source 
of the current layer and by the radiation coming from above, and the interaction 
between the current layer and the layers above are taken into account. There are six 
quantities calculated. They represent the partial upward and downward radiance 
flelds, and incident-from beneath reflection and transmission operator, respectively. 
The precise physical meanings of these quantities can be understood from the follovring 
recursive formulae and the initial values for them. 

Vf - [/ -r(T/»,,T/) [E"‘(t/.T/*|) ^ r(T/MT/) V^-i] 

Vt = Vi-i + REt-x Vf 
vr*E"(T/,T/„)^-f(T/»„T/) V/ 

H, = [/-r(Tj„.T/) REi^xY^ f (t/,t/m) 

Gi = /?£/-, H, 

REj = r(To,T/»i) = r(T/.T/^i) + Gj 

where /sO, ,K-1 and the initial values for them are given by 
Vo = E"(to,t,) + r(T„To) i*^(To) 

Vi = L*(to) 

VS = E*(to,t,) ■» t(T|,To) X.*(to) 

REo = r(To.T|) 

C7o = 0 

Bo = f(To.Ti) 

Z. Calculate the radiance at the atmosphere - earth interface. 

From the forward pass we have calculated the total contribution of the internal 
sources and Uie radiation coming from the top boundary to the downward radiance at 
the surface level, which is . Using the interaction principle again, we can come up 
with the final result for the downward radiance at the earth surface. 

/■•(T*) ♦ 

REk-, 

The upward radiance at this level can be obtained from the lower boundary condition. 

3. Compute recursively the downward and upward radiances from bottom to top. 
This is so called backward pass. 
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L*{ri) = K/ + Q 

L‘(t/) = Vf 4- Hi Z.~(T/ti) 

If only the upward radiance at the top of the atmosphere is needed, an cdtemative 
procedure can oe Lcken. First, calculate r(z,s), r(s,z), t(z,s), f(s.z), and £*(z,t) 
successively for the combined layer, each time adding one mure layer. Finally, obtain 
the corresponding quantities for the entire atmosphere. 

Then, calculate the upward radiance using 

^■(to) = [»-(Tjf,To) «(To.Tjr) Ko (/ -r(To.Tjr) Rc)~^ Mtit.To)] i* Tq) + 

I^'(To.Tir) + f (To.Tjr) Ro U “’•(to.t*') rc)-‘ x 

[S"(ro.Tir) + r(To,TA-) [tB{Tc) ^ ^e /rOic)]] 

2.S. Cooimenta 

Ti;ere are several advantages with this model. 

1. It considers all aspects of the radiation - medium interMtlon, including multi- 
ple scattering occurring within the atmosphere and between the atmt^phere and the 
underlying medium. The complete consideration of all the aspects in multiple scatter- 
ing is especially important when the surface has a high reflectance. 

Z. For surfaces with moderate area size, the algorithm can be used *o model the 
radiance field above it. The modeling accuracy will be good as long as the plane paral- 
lel conflguratlon is valid for the situation. 

3. This algorithm can also be used to calculate the averaged aUnospheric correc- 
tion amount. When this amount is removed from the remotely sensed data, the result- 
ing radiance values will be good indication of the local surface variation. 

4. The algorithm wil! also be useful for removing cloud interference, as long as the 
clouds are not opaque in the wavelength concerned. 

5. Multiple sources, including thermal source, direct solar beam source, and pos- 
sibly even specular pseudo source, are taken into account. 

6 If the upward radiances at the top of the atmosphere are measured systemati- 
cally at a set of vlewii g angles, then the angular radiation fleld at the surface level can 
be inferred by inversion. 

The disadvantages of this r.igorithm include the relatively high computation time 
consumption and relatively larf^e number of input parameters. 

a SOLVINGTHEATVOfiPHEiav:CORRBCT10NPROBl£IIFDRANlNHOI10GEmSQU3SUR- 
FMJB USmC A MONTE CAIOO iinMOD 

a 1. IntroduetioD 

In remote sensing, Landsat satellites look at the earth at near nadir position and 
with a very narrow range of viewing angles. Therefore, by fair approximation for our 
present atmospheric correction algorithm, Landsat looks straight downward at the 
earth surface. 

According to the interaction principle [Orard and Hunt. 1969], the upwelling radi- 
ance for a layered structure atmosphere at the height of the satellite can be 
expressed. 

L'iro) - r{T,,Tc) L*{tc) 4 f (to.t,) L'(t,) 

To=0 is the optical depth at the satellite location, and Tg is the optical depth at the 


ORIGINAL PAGE IS 
- 14 - OF POOR QUALITY 


ground level. The layer concerned is bounded by the levels at Tq. and Tg. r(T^.To) is 
the reflection matrix of the atmosphere of the layer of atmosphere when the incident 
light L*{to) comes from the top of the layer, f (to.t^) is the transmission matrix when 
the incident light L~{Tg) comes from the bottom of the layer, and 2 "(to.t,) is the con- 
tribution of the source within the layer to the i^welling radiance at the satellite level 
To- 

Note that L*(to). L~{Tg), and E~(to.Tj) are vectors describing the angular distribu- 
tion of the radiance, while r{Tg,To) and <(to.t^) are square matrices of corresponding 
size. Also, note that the incident radiances at boundaries are independent of location 
by the above formula. 

For a given atmosphere, the terms r(T^.To) A*(to) and 2 ~(to.Tj) can be calculated 
from an atmospheric radiative transfer model accounting for both zenith and azimuth 
dependence in a absorbing and scattering layered atmosphere. 

The element of I~(to) at the nadir viewing angle, denoted as Z.~(To,'dvf^v)> 
obtained by the satellite measurement. The main concern in remote sensing is to get 
the quantity L~{Tg.iiy,^), vdiich is an element of L~{rg) and where iJv and Vv are the 
zenith and azimuUi angles of the viewing direction of the satellite. For an optically thin 
atmosphere, the term L~{rg,‘&^,fpy) is the leading one contributing to the signature 
I”(To.‘dv.|»v)- If the pattern of the surface radiance is known or a Lambertian property 
is assumed, the surface towelling radiance can be derived given the satellite measure- 
ment and the atmospheric profile, using a layered atmospheric radiative transfer 
model. 

However, for most remote sensing, the situation is more complicated, since the 
upwelling radiance at tho earth surface is not independent of location. We can still 
consider the atmosphere as a layered medium, and within each layer homogeneity can 
be assiimed. In this ceise, the terms corresponding to r(T^,To) />^(to) and S“(To,Ty) are 
still independent of location and remain the same as for a homogeneous surface. 
Therefore, these two terms can be calculated from the layered radiative transfer 
model assuming the underlying surface is black and non-emitting. The only problem is 
that the contribution of the surface upwelling radiance to the iq>welling radiance 
detected by the satellite sensor cannot be simply obtained from the term 
t{TQ.Tg) L~{Tg), because L~{Tg) is now location dependent and the corresponding 
transmission function is also different from the term 1 (to,t,) which is related to the 
integration of the contribution of ground to the radiance at the satellite level over an 
infinite homogeneous surface area. 

Because we are dealing with 2-dimensional surface, we can express the satellite 
measurement by the following equation: 

(To.'^Jvf^v ~ -Artm(To>^vi5®i;) A-L (Toil^vf^wf*^) 

+ L-{Tg,‘&^.ip^.x.y) 

= cosi>v and z ,y are the coordinates for the ground position for the measured pixel, 
aid (.rj are the horizontal coordinates at the satellite level. The first term at the right 
hi.nd side is the element of Z«m(To) at viewing angle and the latter is obtained 

frc m the layered structure atmospheric radiative transfer model. 

Aiim(To) = t(t,,To) Z*(To) + 2"(To.T,) 

If we can estimate then we can obtain L~{Tg,%,^v>(>v) using 

remote sensing measurement and the layered atmospheric model. The Monte Carlo 
method is one of the Outtable candidates for this purpose. 
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3.2. Monte Caiio ^^roadi 

Only the transmission problem is treated, because the reflection and source prob- 
lem can be solved by the previous model. For a s}rstem composed of a layered atmo- 
sphere and an underlying inhomogeneous surface, a combination of the azimuth depen- 
dent radiative transfer model for the L^tm term and the Monte Carlo method for 
transmission of surface upwelling radiance through the atmosphere layer cam give a 
complete expression of the satellite measured radiance. 

The proposed model is based on the fact that the contribution of the surface 
upwelling radiance to the signature of a certain pixel is made by not only the ground 
element corresponding to the pixel concerned, but also the surrounding area as well. 
By the reciprocity principle, the pattern of the contribution made by the ground area 
can be mimicked by a reverse process, a process in which a beam of photons impinges 
at a given point at the top of the layer and Anally some of them get to the surface and 
make a spread pat .ern the so called "point-spread" function. In other words, if we can 
And the point-spread function for a beam impinging on a given atmosphere as a func- 
tion of incident angle, we can obtain the quantitative expression for the pattern of the 
contribution of surface upwelling radiance to the satellite measurement. This relation 
can be written as the following formula; 

f (To,f.T?,'dv.l»v: r,.x.y.4,v) - f(T,,x,y,‘d,<p:To,t»?,dv,<Pv) 


The is the basis of the present algorithm, where t is 

_ AI(T|,Xi,y„d„5P,:T2,X2,y2.'l>2.9i>2) 

fiT„X,,y,,1>„^.:T2,X2,y2,d2,«»2) - 2(T2,X2,y2.^.«»2)COSt>2d«8 0M2 


tL is the increment of radiance at location (T|,x,,y,) and in the direction con- 

tributed by the radiant flux at location (T 2 ,Xg,yz) in the direction 

Note that for this process, the conservation of photons does not hold because 
some of the incident photons may be absorbed and some of them may get out of the 
layer through the uppi boundary. In either case, the photons are lost and have noth- 
ing to do with the transi .ssioi. point-spread function. Besides, in this configuration, we 
need not consider the emission problem because the previous adding model has taken 
care of that process already, and this makes our Monte Carlo method relatively simple. 

We need to determine the point-spread function using the Monte Carlo method. 
Several conditions are given: 

1. The beam of photon impinges at the top of the layer of atmosphere concerned. 

2. The atmosphere has layered structure and the profiles, including molecular 
absorption emd scattering, Mie scattering components, and temperature, are known. 

3. For each sublayer, homogeneity is assumed with the scattering and absorption 
properties assumed to be the average value derived from the given profile. 

4. The relation between optical depth and height should be derived from the 
atmospheric profile, and be used throughout this algorithm, since this knowledge is 
crucial for calculating the horizontal position at which the photon gets out of the atmo- 
sphere. 

The essence of the Monte Carlo method is that the scattering and absorption of the 
photon can be statistically simulated by a sequence of photons randomly colliding with 
the mediuiT particles before they Anally get absorbed or escape from the medium con- 
cerned. Ef.ch scattering or absorption is a random event of collision. However, the 
general trend is governed by some statistical rule according to some probability func- 
tion of the processes. 

In order to mimic this random process, we need some quantities based on some 
statistical rule. These are the distance traveled by the photon between two random 
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collision events, the direction of each traveling path, and the chance for the photon to 
be absorbed iq>on collision. Also, we need to keep track of position of the photon in 
three-dimensional space in order to know whether the photon is in the atmosphere, in 
which .sublayer the photon is, and how far it is from the original horizontal position. 

The following sections give the mathematical expressions of these events and 
quantities. Similar descriptions can be found in some representative papers [Chs/itueU 
and Everett, 1959 ; House and Avery, 1969 ; Pearce, 1977 ]. 

3.2.1. Free Path Length 

The distance traveled by a photon between collisions is called the free path length. 
When it is measured in the same manner as the optical depth, it is dimensionlei'is and is 
called optical distance. The probability density function of the occurrence of noncolli- 
sion for an optical distance I is 

P(0 = e"‘ 

The probability that no collision occurs in the range of optical distamce from 0 to I is 

i 

r, = f p{l') dt' = 1 - e~‘ 

0 

This equation sets up a unique relation between a random number t\ (O.O^ri^l.O) and 
an optical distance I , and I can be solved explicitly as 

I = -ln(l -r,) = |ln(l -r,)| 

Iheretore, the distance between scatterings for each path of each photon can be deter- 
mined according to a random number T] generated in some manner. 

3.2.2. Chance of Absorption 

The single scattering albedo S of the current sublayer gives the probability 1-0 
that the photon is absorbed upon collision Conceptually, we can use the idea of photon 
bundle instead of single photon. Suppose initially a photon bundle impinges into the 
atmospheric slab. After collision, only a fraction of the photon bundle keeps going and 
a portion is absorbed. For each collision, the current fraction is multiplied with the 
single scattering albedo of the current sublayer to get the fraction of the photon bun- 
dle continuing. Such a process is repeated until the remaining portion of the photon 
gets out of the slab or is less than some threshold. 

3.2.3. Direction of Scattering 

For each scattering event, two independent emgular variables can be obtained 
from random processes, the scattering angle 9 and the azimuth angle $ which is meas- 
ured in the plane perpendicular to the original travel direction 

The scattering angle 9 is determined in the following way: 
define 

e 

F’(9) = f p{9')siD9'd9' 

0 

9' is a dummy veu*iable. p{9') is the phase function for scattering angle 9’ for the 
current scattering sublayer. Because the integration of p(9') over the range from 0 to 
n is 2, we need to multiply the quantity expressed by a factor of 0.5 to normedize it. 
Then we can relate such a quantity to a random number rz generated for the reuidom 
scattering process to determine the scattering angle 9. 

= re 
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The angle $ is within the range from 0 to 2rr, therefore 
♦ = 2tt rs 

where rg is another random number. Both rg and rg are within the range from 0 to 1. 

Knowing the original travel direction and the scattering angle and azimuth 9 
and 4, the travel direction for the next path, be determined from 


cosdg = cosi9) cos9 - sin'di sin9 cost 

_ cosd|Siny|Sin9cost » cosyiSinSsint sintfiSin^giCosB 


tan(0g = 


cosiliCos^iSin9cos# - sinf^isinOsint -i- sini>iCOS{»iCos9 


These relations can be derived using the formulae for 3-dimensional axis rotation. 


3.2.4. Distance of PenetraHon of Photon in the Slab 

In terms of optical depth, the penetration is determined by 
At = 1 cosdg 
Tg = Ti + At 

where I is the optical distance, T| and Tg are the optical depths fo^ two successive colli- 
sions, and At is the increment of optical depth between the two collisions. 

The height at which the collision occurs can be calculated from the relation 
between the optical depth and the height according to atmospheric profile. Suppose 
the heights for two successive collisions are hi and hg, the distance traveled between 
collisions is 

d = (hg - hi) / cosdg 


3.2.5. Horizontal Displacement 

The horizontal displacement can then be calculated. 
bz = d sindg COS^g 
Ay = d sindg sin^Pg 

where A* and Ay are the coordinate increments in S-N and E-W directions, respec- 
tively, and the horizontal location xz.yz can be derived given the location before the 
current travel pathxj.yi. 

Xg = 2i + Az 

yg = y, + Ay 


3.2.6. Escape 

If the current accumulative optical depth Tg is less than 0, the photon is reflected, 
and this case can be ignored since the present algorithm only deals with the transmis- 
sion problem. If Tg is greater than the total optical depth of the slab, the photon 
escapes at the bottom of the slab and contributes to the point-spread function. The 
horizontal location should be calculated using a modified traveling distance which is 
shorter than what has been calculated from the random process. The horizontal loca- 
tion thus obtained and the current travel direction should be recorded properly in 
order to get an accurate uesciiption of the point-spread function. 

This concept can also be slightly modified using the photon bundle mention before. 
Suppose a fraction of a photon bundle fi travels in direction starting from T,x,y. 
For d>rr/2, the fraction of the photon bundle /g that will escape at the top of the slab 
is 
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/j = e-(T/|co.*|) 

For 0<i><jr/2, 

f2 = f\ 

For i>=tt/2. 
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/z = 0 

T* is the total optical depth. The horizontal position of escaping in the manner 
described earlier. The remaining portion of the photon bundle / j -/ g will still travel 
within the slab. The distance traveled related to random number r, can be determined 
from 


I = |ln(l - r, (1 -^^)| 

J I 

Such a process along with the absorption process upon collision for non-conservative 
scattering makes the remaining portion of the photon bundle smaller and smaller and 
eventually the remaining fraction is so small that it can be neglected. Such a treat- 
ment can save computation time for given precision. 


3.2.7. IXscretizaHon of Ar^lss and Distance 

When the photon travels within the slab, the locations of the photon are calculated 
in a continuous manner so that no error is accumulated. However, at the point the 
photon escapes at the bottom of the slab, we have to record the escaping location and 
direction of the photon in a discrete manner. Four dimensions of discretization may 
be applied since the distribution of the escaping photon is a function of location coordi- 
nates of x.y and angular coordinates An array of 4 dimensions takes a lot of com- 
puter memory space, and to make the recorded photon number in each element of 
such an array statistically meaningful, the total number of photons neeaed in the 
experiment must be tremendously large. Therefore such an algorithm is infeasible. 
Under the assumption that the photon impinges at the top of the slab perpendicularly, 
the point-spread function is axisymmetric and we need to record the information along 
a single radius only instead of over the entire x.y plane. Therefore a recording array 
of only 3 dimensions is needed. We fix the single recording radius along the x axis emd 
put on it the results from the location away from it by rotation. The rotation angle 

IS 


tan^p, = yz/i2 

Note that the signs of Xg and yg allow <pr to be determined over the range 0 to 2n. After 
rotation, the zenith angle remains the same while the azimuth angle is changed. 
Therefore the new zenith angle d is 

d = dg 

and the new azimuth angle ^ is 
<P = <PZ-Vr 

The distance of the total horizontal displacement is 
r = Vxg* + yg* 

Later on we will see that for a Lambertian surface, the discretization of d and ^ are not 
necessary, only the spatial discretization is needed. For such a case, the computer 
resources needed can be reduced, and the non-axisymmetric case would be computa- 
tionally possible. 
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For the axisymmetric case, the discretizations of r , i), and ^ can be performed in 
the following manner: 

For the r-dimension, an equal increment discretization scheme is carried out. 
The increment of r should be equal to or less than the ground resolution of the pixel. 
For zenith angle ‘d-dimension, the double Gaussian quadrature rule for cosine of zenith 
angle is the suitable one because the result obtedned will be comparable to those pro- 
duced in the layer adding radiative transfer model. For the azimuth angle domain, the 
best discretization is the equal increment scheme. The range covered by each element 
can be determined by finding the midpoints between the adjacent discrete points. 


3.2.B. Calculation of the Point Spread Function 

In digital image processing, we need to express the point-spread function {PSF) in 
matrix format. If we superimpose such a matrix over a matrix whose elements 
represent surface irradiance, then the product of two corresponding elements in these 
two matrices represents the contribution made by the irradiance from that particular 
ground element to the radiance measured by the satellite, and the sum of the products 
within the PSF window is the total contribution made by ground radiation to the signa- 
ture of the pixel that coincides with the central element of the PSF matrix. 

In order to construct such a matrix, we need some corrections to the results 
obtained directly from the Monte Carlo me^od. 

First is the correction for the difference between areas covered for a fixed incre- 
ment of radius Ar but for different values of radius r. This is obvious, since the central 
point's influence ranges from 0 to }>^r and therefore the infiuence area is (n/4) Ar^. 
Other points’ infiuences range from r-J^r to r+}^Ar with area 2rrr Ar. Therefore it is 
necessary to make the area correction to get the contribution per unit area. 

Secondly, the difference between the area of a circle and the area of square or 
rectangle must be corrected. For a matrix, the most convenient shape for each ele- 
ment is square or rectangle. 

After these corrections, the contribution made by the surface radiation from a 
unit area of Lambertian surface for the central ground element and the surrounding 
elements can be expressed as 


tuo = * * 


E E Nirc.A<p) C 


Nf 


E E C 


= -i- 


M ^ 


Uq is the central element of the matrix of the point-spread function, and ttv is an ele- 
ment other than the central one with distance from center r. N is the 3-dimensional 
array obtained in the Monte Carlo algorithm, and Nt in the total number of photons 
incident at the top of the atmospheric slab in the experiment. C corrects for the 
difference between the area of a circle and that of a square or rectangle. For a square, 
C = 4/ rr. From these equations, it is obvious that only the radius dimension is neces- 
sary for Lambertian surface. 

The convolution of such a point-spread function matrix with the matrix of ground 
radiance, in which each ground element is Lambertian, gives the total contribution of 
surface upwelling radiation field to the satellite measurement, and the satellite meas- 
ured radiance can be expressed as 

= LHniTo.^v><Pv) 
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+ 22 A(x*H.y,_^)L*(T,.tfv.|Pv.*(.V/) 

i>0i>0 

h is the point-spread function matrix of size A/xAf ; It is convolved with a ground matrix 
of size NxN. The elements of h outside the AlxA# window size are all zero. The last 
term on the right hand side is the convolution term. 

The point-spread function matrix thus defined only works for the isotropic surface. 
For anisotropic surface, we need to define the contribution factors for the central 
ground element and the surrounding elements, respectively. 

~ ** A#(t,.xo.Vo) Nt 


= S ETA(T„.i^.<fi.Xo.yo) fV(ro.t>.^) 

E E C br 

t _ » e 

^ B Af(Tj,X(,yy) Nt 


ETA is the anisotropic factor for upwelling radiation under given illumination condition, 
and is defined as 


£Ti 4 (T,.d,^,x.y) 


(Tg,d,y,x,y) 

M{Tg,x.y) 


M is the exitance at given location. Then the satellite measurement of the radiance is 


+ 22 M{Tg,Xi,y^) 

i«0 isO 

/14 is the element of the matrix composed of the contribution i 4 ^'s. For a Lambertian 
surface, a simpler formula can be used. 


* t t MX. -i.yi-i) L~{Tg,Xi.yi) 

iso j*o 


h is the element of the matrix composed of weight wq and vv's. 

From these equations, it is obvious that if the point-spread function is derived 
from the Monte Carlo method and if the surface exitance and the pattern of the radi- 
ance angular distribution at each surface location are known, the vertically upwelling 
radiance at satellite level can be determined. 


This is the forward aspect of our problem. However, in remote sensing, the 
reverse aspect of the problem is of more significance, but it is more difficult, too. In 
the following section, we discuss the inverse problem. 
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3.3. Retrieval of the Surface BliUaoe 

After removing the atmospheric reflection and source term from the satellite 
measurement, we have an image g. Also, we can construct a matrix H from the point- 
spread function such that the size of H is comparable to the lexicographic form, by 
which the rows of a matrix are stacked to form a vector of the matrbc g , and the lexico- 
graphic form of the corresponding matrix of surface exitance /. The relationship 
between f ,g, and H is 

g^fif 

If there is no noise other than the atmospheric effect, then the surface exitance can be 
retrieved by the inverse Alter given g . 

The Fourier transform technique can reduce the computation time for the inverse 
Altering [>bidreii« and Hunt, 1977], In the presence of noise, other types of Alters may 
be employed using the point-spread function along with some a priori knowledge of 
noise [Andrews and Hunt, 1977]. However, those procedures are more complicated 
and the extent to which the result can be improved strongly depends on the 
signal /noise ratio. The Fourier transform technique similar to that used in Inverse 
Altering c an also used in those Alters, to reduce computation time. 
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